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Abstract —Current and future radio interferometrie arrays such as 
LOFAR and SKA are characterized by a paradox. Their large number 
of receptors (up to millions) allow theoretically unprecedented high 
imaging resolution. In the same time, the ultra massive amounts of 
samples makes the data transfer and computational loads (correlation 
and calibration) order of magnitudes too high to allow any currently 
existing image reconstruction algorithm to achieve, or even approach, the 
theoretical resolution. We investigate here decentralized and distributed 
image reconstruction strategies which select, transfer and process only 
a fraction of the total data. The loss in MSE incurred by the proposed 
approach is evaluated theoretically and numerically on simple test cases. 

I. Introduction 

Since the commissioning of the first large radio interferometers in 
the 70s and 80s (such as the VLA in the USA and the WSRT) radio 
astronomy in the range of large wavelengths has grown dramatically, 
particularly with the development of more and more extended antenna 
arrays. In the prospect of the most sensitive radio telescope ever 
built, the SKA which will he operational in the 2020s, several new 
generation radio telescopes are being built or planned (LOFAR in the 
Netherlands, ASKAP and the Murchison Widefield Array Australia, 
e-MERLIN in the UK, e-EVN based in Europe, MeerKAT in South 
Africa, JVLA the United States). 

As an example, LOFAR consists of 48 groups of antennas (sta¬ 
tions), among which approximately 35,000 elementary antennas are 
located in the Netherlands. The “superterp”, the heart of LOFAR 
is a super-station: a cluster of six stations. Eight other stations, 
totalizing approximately 13,000 antennas are located in the sur¬ 
rounding countries. A project of a new super-station in Nanjay 
(F) is under consideration. Within each station, antennas form a 
phased array which allows for digital heamforming simultaneously 
in several directions and frequency bands. The heam-formed data 
from the stations are centralized at the University of Groningen 
in the Netherlands where a supercomputer is responsible for the 
combination of the heam data from all stations. The resulting data 
are then stored on a cluster of ASTRON, the Netherlands Institute 
for Radio Astronomy, where the images (and other deliverables) 
are reconstructed. As a mean of comparison SKA will totalize 2.5 
millions antennas, with a square kilometer collecting area distributed 
over an area of « 5,000 km diameter. 

Beyond specific objectives that distinguish these new fully digital 
“software telescopes”, they are all characterized by a great flexibility. 
Another common point is the amount of data which must be trans¬ 
ferred to the central computer and processed. It amounts to 1 ter¬ 
abit/second for LOFAR and will be of the order of 14 exabyte/day 
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for SKA (more than 100 times the global internet traffic). LOFAR 
uses a 1.5 Blue Gene/P for the data reduction and the computation of 
correlations. IBM et ASTRON will develop by 2024 a supercomputer 
to process and store 1 petabytes of data everyday (2). 

This correspondence investigates the possibility to distribute the 
image reconstruction over the super-stations. The main objective is 
to avoid centralization of the sampled electromagnetic fields acquired 
by all stations in order to reduce the data transfer and the exponential 
increase in the calibration and computational load. 

Section II recalls the basis of radio astronomy with aperture 
synthesis and proposes a strategy where each super-station uses all 
its antenna and one reference signal from the other super-stations. 
The loss of performances that follows is evaluated on a simple 
model using the Cramer Rao Lower Bound (CRLB). Section III 
shows that the image reconstruction problem can be written as a 
global variable consensus problem with regularization. Numerical 
simulations illustrate the performances of the proposed approach. A 
concluding section presents perpectives. 

11. Aperture synthesis for radio astronomy 
A. Standard aperture synthesis model 

This section provides the basic equations of radio astronomy 
with multiple sensor array and describes a partial aperture synthesis 
strategy which aims to reduce data transfer, allowing a decentralized 
image reconstruction. 

To simplify the notations and without loss of generality, we will 
not make explicit the wavelengths dependence and the Earth rotation 
and assume punctual antennas. The coordinates of the stations (within 
each station, a beam is created from the phased array) in a plane 
perpendicular to the line of sight are denoted as Vj and the map 
of interest (the “image” of a region of the sky) is x{p) where p 
denotes the angular coordinates on the sky. The fundamental equation 
of interferometry relates the Fourier transform of the map to the 
spatial coherency (visibility) of the incoming electromagnetic field. 
A measurement of the coherency is obtained by correlating the signal 
acquired by a pair of stations {i,j) properly delayed located at Vi 
and Tj, giving in the noiseless case a point of visibility at spatial 
frequency ui = rj — rp. 

v{ui) = J x{p)e~^^^'^^^dp (1) 

See for example 1101 for a comprehensive description of radio 
astronomy and signal processing related tools. 

Computation of v(ue) obviously requires the transfer of signals 
from stations i and j in the same place. The stations are nor¬ 
mally grouped in “super-stations” (e.g. the superterp for LOFAR) 
accounting for low frequencies (||ri — small). Resolution is 


then increased by correlating signals between stations that can be 
located up to thousands of kilometers from each other. 

B. Reduced synthesis aperture model 

In general, all possible correlations are computed in order to 
maximize the (u, v) coverage (set of uj). This requires a centralized 
system architecture. The resulting visibilities measurements v{ue) 
are associated to the filtering of the visibility function v{u) by the 
Full Aperture (FA) spatial transfer function 

Af{u) = S{r - rfc) * ^ ^(r - rz)^ (w) (2) 

M 

= ^ 'Wi8{u — u^) — Al{u) + Ah{u) (3) 

where * denotes the convolution, the weights we count the number 
of beam pairs measuring the same spatial frequency and M is 
the number of different sampled frequencies. The term Al(u) is 
associated to “intra-super-stations” low-frequency correlations and 
A^(u) to high-frequency “inter-super-stations” correlations. Note 
that in order to simplify the derivations, we will also denote super¬ 
station a single remote station. 

In order to reduce the amount of transferred data, we propose to 
investigate a solution which consists in: 

1) exploiting all the low-frequencies by computing locally all 
the correlations inside each super-station. In this case, the 
low frequency term associated to the spatial transfer function 
(denoted as Ar{u)) is still Al{u). Denote as Sk the set of 
indices associated to the beams in super-station k. 

Al,l{u) ^ ^ 5{u - {Vm - rn)) + S{u + {Vm - Vn)) 

m<nGSi 

al{u) = y^^Aijju) 

I 

2) Recovering the high frequency information by transferring only 
one single beam signal from each super-station to all other 
remote super-stations. If Ck is the index of the reference beam 
in super-station k transferred to the other super-stations, the 
resulting sampling pattern of visibilities associated to super¬ 
stations k and I is, see Fig. 

Ak,i{u) = ^ 5{u - (rm - rck)) 

m^Si 

-^h{u) = y^^Akdju) 

k^l 

This solution is motivated by the fact that it allows to correctly sample 
the high frequencies, at the cost of a reduced SNR, while reducing 
the number of transfered electromagnetic signals. Other strategies 
than transferring a single antenna signal are of course possible. A 
solution which aims to preserve the SNR is to replace the signal 
indexed by Ck by averaging of neighboring beams. It is important to 
emphasize that unless the averaging process is taken into account in 
the measurement equation, the bias introduced by the averaging of 
nearby beams must be negligible. This gain in SNR will be denoted 
as p > 1 in the sequel. Finally, the overall Reduced Aperture (RA) 
antenna spatial transfer function is: 

Ar{u) = Al(u) + pAh(u) (4) 

This strategy reduces the transfer of beams data w.r.t. a centralized 
processing as long as the number of stations inside each super-station 



Fig. 1. Reduced Aperture (RA) synthesis using super-stations k and 1. 

is larger than the number of super-stations: for Nsa super-stations of 
Ns stations each, the first requires Nss{Nss — 1) transfers whereas 
the second NssNs- 

C. Analysis of performances on a simple model 

In order to evaluate analytically the loss related to the use of RA 
synthesis w.r.t. to a FA we consider a simple one dimensional case 
where the map is a shifted Gaussian shape with flux a: 

Q, _(P-PO)^ 

*(p) = —7=e (5) 

TjyZTV 

The unknown parameters are 6 = (a, 77 ,po)- The visibilities are: 

v ( ui ) = + ne , £= 1 ...M (6) 

where ne is a measurement noise assumed independent Gaussian 
circular with ni ~ Afc{0,wf^a^). The coefficient takes into 
account the variance reduction that occurs when the visibility v{ui) 
is estimated from we different baselines. 

The elements of the Fisher information matrix /(0) are computed 
using the Slepian-Bangs formula p. 293] which gives: 

2 f 

1(9) = — -2FaS2 
\ 0 
M 

1=1 

Note that I{6) is not a function of the source position po. 

We compare the CRLBs on a, rj and po for two spatial transfer 
functions. In both cases the aperture configuration consists of two 
uniform sub-apertures separated by D in order to sketch the behaviour 
of two super-stations: 

rk = -D/2-kA, k^0...L (9) 

VL+k+i = D/2 -I- kA, k ^ 0.. .L (10) 

where D > {L + 1)A. In the FA mode, and for u > 0: 

L 

Al{u) = ^ 2(L + 1 - £)5{u - £A) 

L 

AUu) = ^ (i + 1 - \£\) S{U -{D + {£ + L)A)) 

e=-L 

L is assumed even, L = 2q and we consider in the RA mode 
that Ck = ±q: the reference beam is in the middle of the opposite 
super-station. As noted above, the low frequency term Al{u) does 
not change. The high frequency spatial transfer function is now for 
u > 0: 

L 

Al{u)^Y.^{u-{D + {£ + q)A)) 

1=0 


-2FaS2 0 \ 

0 (7) 

0 'iFa^S2 ) 

( 8 ) 





FA synthesis spatial transfer function: Af(u) 



where the regularization term ^(a:) can be for example a sparse 
analysis El, a sparse synthesis lITTl or a hybrid prior (4). 

Denoting T^'v the zero-padded visibilities and y = F^T^v the 
so called “dirty image”, (Uj can be rewritten as 


y = Hx + e 


(13) 


where H — F'T'WTF is a (circulant) convolution matrix. Using 
yz ||t1z||2 = ||: 2:||2 and F"^F — I, we have ||i; — Ga :||2 = \\y — 
Hx \\2 and the inverse problem il2i turns to be equivalent to the 
deconvolution problem 


Fig. 2. ID illustration of Af{u) and Ar(u). Bottom plot shows |u(u)| for 
2 characteristic values of rj. 



Fig. 3. CRLB for the FA: CRLBjr(-), and the RA: CRLBfl(-). 

The source width rj plays a central role in the estimation. For a 
point source, y —>■ 0, high frequency measurements will bring a 
lot of information while performances for a very extended source 
(rj —>■ oo) will be independent of the inter-stations visibilities. Figs.|^ 
and give the results obtained for a configuration defined by L = 4, 
D = 10 and A = 0.1. The source parameters are a = 1, po = 0 
and results are given for different values of rj. The gain p is fixed to 
p — 2. Fig. shows the two spatial transfer functions Af{u) and 
Ar{u). Fig. 1^ shows the CRLBs associated to Af{u) and Ar{u), 
denoted respectively as CRLBf and CRLB/j. The thresholding effect 
when rj « 0.05 reflects the shape of the visibility modulus given 
in Fig. 1^ the frequency contribution of the source at the inter¬ 
station baseline becomes negligible for y > 0.05. For y < 0.05 
the order of magnitude of the loss of performances is 4dB. This loss 
of performance naturally strongly depends on p, e.g. when p — 1 the 
loss is 13dB. 

III. Distributed image reconstruction with partial 

APERTURE SYNTHESIS MODEL 
A. Decentralized map reconstruction 

The classical model for the map reconstruction is obtained vec¬ 
torizing the sampled map x £ (R^)^ and the visibilities v £ 

M < N, and reads El, El: 

v^Gx + n,G^WTF (11) 

where F is the N x N Fourier transform matrix, W is a diagonal 
weighting-matrix including various operations (calibration, signal to 
noise weighting), and T is a 0/1 binary M x N matrix which codes 
the sampling of the visibilities in the frequency plane. The noise 
vector is assumed n ~ With these assumptions, the 

image restoration problem is usually casted as the general inverse 
problem: 

a: = arg min \\v — Gx\\ 2 -\-Q.(x) (12) 

a,g(l.+ )N 


X = arg min lly — iTailli + D(a;) (14) 

a:6(K+)™ 

The vector v can be partitioned as: 

,...), vl. = (itl,i,ul,2,---) (15) 

where denotes the vector of visibilities obtained by correlating 
all the beams inside super-station k and vector Ufe,;, k ^ I, hy 
correlating beams in station k with reference beam indexed by ci 
in station 1. More specifically Vk,k is associated to Ak,k(u) and Vk,i 
is associated to Ak,i{u), see sec. II.B. Using the same partitioning 
for G, W and T, Eq can be rewritten as : 

A = arg min llufc . — Gfea;||i-F D(a;) (16) 

a:g{R+)JV ^ 

with Gk = WkTk- Using again the property of the zero padding 
matrix tJ, inside each norm in the sum |l6|l, we obtain: 

a; = arg min \\yf. . — Hkx\\l + Q{x) (17) 

a,6(E-l-)NZ^ 

where Hk = F^T^WuTkF. Note that H = Y^^.Hk and y = 
Y^kVk- 

In l |16^ , each sub-problem in the sum amounts to reconstruct x 
from the intra-super-station and inter-super-station visibilities Vk,-- 
In (|17|, j/j, , is a dirty image obtained using only the visibilities Ufc,. 
and each sub-problem amounts to reconstruct x from .. 

Eqs. jl6|17^ are particularly interesting for the derivation of 
distributed optimization algorithm, since they correspond to a global 
variable consensus problem 0, Q. The next subsections evaluate 
the impact of the partial aperture models w.r.t. the standard model by 
numerical simulations. 

B. Array configuration and aperture synthesis 

The shape of the sensor array used in the simulations consists 
of 10 super-stations of 10 stations each, as shown in Fig. The 
stations follow a classical Y configuration with different rotations. 
Ten measurements are performed in a range of 3 hours taking into 
account the Earth rotation. The RA is obtained by computing for 
each station the correlation with the center beam of the other super¬ 
stations and with p — 2. For this experiment the ratio D/{LA) is 
20 and the visibilities are binned in a 256 x 256 grid with Shannon 
sampling. The corresponding (u, v) coverage is shown in Fig. 

The original image is shown in the upper part of Fig. along 
with its Fourier transform. Note that the image contains both low 
and high frequency. The noise of variance is such that the 
measured visibilities SNR is 70 dB. The observed dirty images are 
also reported. The RA clearly leads to a less detailed dirty image. 
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Fig. 4. Left: sensor array configuration (without Earth rotation). Final (u, v) 
coverage in normalized frequencies for FA, Eq. ^ and RA, Eq. |^. 


Original image Original FT (log scale) 



Fig. 5. Original image and its Fourier transform (up). Dirty images observed 
with the full aperture (down left) and reduced aperture (down right). 


C. Distributed image reconstruction 

The image reconstruction is performed by solving problem l |17^ 
with a regularized global variable consensus ADMM algorithm as 
described in m. The regularization term is Q,{x) = A||a ;||2 with 
A = 10“®. This regularization ensures that the problem is strictly 
convex and limits the bias. At each iteration, this algorithm requires 
to solve a large scale linear problem of size the number of pixels 
N in each super-station, [T] Eq. (7.6)]. This linear problem can be 
easily solved in the Fourier domain as discussed in O. Note that the 
quadratic regularization can be included in this step. The consensus 
step (H Eq. (7.7)] is simply a projection on the positive orthant. 
A super-station sends the current image, Lagrangian multiplier and 
receive the consensus image. 

Fig. § shows the reconstructed images obtained by solving 0 
with and without positivity constraints for both aperture cases, along 
with the relative norm of the error e. The FA leads obviously to 
a better reconstruction in both cases. However, while the loss in 
performance is rather important in the unconstrained case (e = 15% 
error instead of e = 23%), this gap is significantly reduced with the 
positivity constraint (e = 9.5% error instead of e = 12%). 

IV. Conclusion and perspectives 

This correspondence investigates a distributed strategy for the 
image reconstruction problem in radio astronomy when the number 
of stations inside each super-station is larger than the number of 
super-stations. It relies on a reduced aperture synthesis where each 
super-station uses all its beam and a single reference beam from 
the other super-stations. Part of the missing information for each 
super-station is then exchanged during the consensus step of the 
distributed algorithm. The loss of performances, related to the use of 
a reduced aperture synthesis, is evaluated on the image reconstruction 
by computer simulations. 



Fig. 6. Reconstructed images and error e. Up (down): reconstruction without 
(with) positivity constraint. Left : FA mode, right : RA mode. 


The approach proposed in the paper processes all the data after they 
have been acquired. A natural extension is to reconstruct sequentially 
the image. Among the benefits of this setup which perfectly fits 
the operating mode of interferometers which progressively fill the 
frequency plane using the Earth rotation, is the possibility to optimize 
in real-time the observation mode. A straightforward solution is to 
make a number of iterations of the reconstruction algorithm of section 
|III-C| after each measurement and using a “warm start”. A much 
more challenging perspective it to select sequentially the frequency 
measurements used at each iteration. Whereas this communication 
relies on a “deterministic” pattern of frequency measurements, a 
better strategy would be to select at each iteration, among all 
visibilities, the ones that optimize the reconstruction according to 
some predefined criterion. 
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